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ABSTRACT 

We have investigated the relaxation of a hydrostatic hot plasma column con- 
taining toroidal magnetic field by the Current-Driven (CD) kink instability as 
a model of pulsar wind nebulae. In our simulations the CD kink instability is 
excited by a small initial velocity perturbation and develops turbulent structure 
inside the hot plasma column. We demonstrate that, as envisioned by Begelman, 
the hoop stress declines and the initial gas pressure excess near the axis decreases. 
The magnetization parameter cr, the ratio of the Poynting to the kinetic energy 
fiux, declines from an initial value of 0.3 to about 0.01 when the CD kink instabil- 
ity saturates. Our simulations demonstrate that axisymmetric models strongly 
overestimate the elongation of the pulsar wind nebulae. Therefore, the previous 
requirement for an extremely low pulsar wind magnetization can be abandoned. 
The observed structure of the pulsar wind nebulae do not contradict the natural 
assumption that the magnetic energy fiux still remains a good fraction of the 
total energy fiux after dissipation of alternating fields. 

Subject headings: instabilities - magnetohydrodynamics (MHD) - methods: nu- 
merical - (stars:) pulsars: general 

1. Introduction 

The pulsar wind nebulae (PWNe) may be considered as a relativistically hot bubble 
continuously pumped by an electron-positron plasma and magnetic field emanating from the 
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pulsar. Pulsars lose their rotational energy predominantly by generating a highly magnetized, 
ultrarelativistic wind. The wind presumably terminates at a strong reverse shock and the 
shocked plasma inflates a bubble within the external medium. The synchrotron and inverse 
Compton radiation from the shocked plasma is observed from the radio to the gamma-ray 
band (see the review by Gaensler & Slane 2006). 

Close to the pulsar the energy is carried mostly by electro-magnetic flelds as Poynting 
flux; however, the common belief is that at the termination shock the wind must already be 
very weakly magnetized. Simple spherical models of PWNe suggest that the magnetization 
parameter a, the ratio of the Poynting to the kinetic energy flux, needs to be as small as 
0.001-0.01 just upstream of the termination shock (Rees & Gunn 1974; Kennel & Coroniti 
1984a,b; Emmering & Chevaher 1987). The reason for the required low value of a at the 
termination shock is that conservation of the magnetic flux in the effectively incompressible 
subsonic flow downstream of the termination shock implies rapid increase in the magnetic 
field strength with radius and the field within the nebula could exceed the equipartition 
value if the magnetization at the termination shock is not extremely small. Extensive ax- 
isymmetric MHD simulations of the fiow produced by the pulsar wind within a plerionic 
nebula (Komissarov & Lyubarsky 2003, 2004; Del Zanna et al. 2004, 2006; Volpi et al. 2008; 
Camus et al. 2009) show that one can account for the morphology of PWNe, including the 
remarkable jet-torus structure, with a ~ 0.01. If the magnetization were larger, the nebula 
would be elongated by the magnetic pinch effect beyond observational limits. Such a low 
value of (7 is puzzling because it is not easy to invent a realistic energy conversion mecha- 
nism to reduce a to the required level. This problem, often referred to as the "cr problem" 
is widely discussed in the hterature (see recent reviews by Arons 2007; Kirk et al. 2009). 

One has to stress that all the available observation limits on cr are obtained from the 
analysis of the plasma flow and radiation beyond the termination shock, where the upstream 
cr is calculated from the ideal MHD jump conditions as if the Poyntning flux is transferred 
by large scale magnetic flelds. However in the pulsar wind, most of the energy is transferred 
by waves, which an obliquely rotating magnetosphere excites near the light cylinder. These 
waves cannot propagate within the nebula because the wavelength (on the order of the light 
cylinder radius) is less than the particle Larmor radii. The above mentioned observational 
limits on a refer only to the mean magnetic fleld remaining after the oscillating part is erased. 

In the equatorial belt of the wind, the sign of the magnetic fleld alternates with the pulsar 
period, forming stripes of opposite magnetic polarity (Michel 1971; Bogovalov 1999); such a 
structure is called a striped wind. In the striped wind, Poynting flux can be converted into 
particle energy flux when the oppositely directed magnetic flelds annihilate. Observations 
of X-ray tori around pulsars (Gaensler & Slane 2006) as well as theoretical modeling of 
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the pulsar wind (Bogovalov 1999; Spitkovsky 2006) suggest that it is in the equatorial belt 
where most of the wind energy is transported. Therefore, in the equatorial belt magnetic 
dissipation of the striped wind is the main energy conversion mechanism in pulsars. It has 
been shown that due to relativistic time dilation, complete dissipation could occur only on 
a scale comparable to or larger than the radius of the termination shock (Lyubarsky & Kirk 
2001; Kirk & Skjaeraasen 2003, Zenitani & Hoshino 2007). However, the alternating fields 
still annihilate at the termination shock (Petri & Lyubarsky 2007). At higher latitudes, 
where the magnetic field does not change sign, fast magnetosonic waves may transport a 
significant amount of energy. These waves can decay relatively easily (Lyubarsky 2003) but 
can release only a fraction of the Poynting flux into the plasma, because at these latitudes 
most of the energy is carried by the mean magnetic fleld. 

The fraction of the energy transferred by the mean field can be found only from 3D 
numerical simulations of the pulsar magnetosphere. Even though this fraction is still not 
known, this fraction is less than 1/2 because the angular distribution of the Poynting flux in 
the pulsar wind is a maximum at the rotational equator, where the mean field is zero. This 
suggests that a becomes less than unity after the waves decay. Therefore, at a quantitative 
level the a problem is partially solved if Poynting fiux is converted into plasma energy via 
dissipation of the oscillating part of the field. However, the residual a still cannot be as 
small as required by axisymmetric models. Therefore the question still remains as to how 
the mean field cr, which is only somewhat less than unity, could become extremely small. 

Since no mechanism had been found for the extraction of energy from a large scale, 
axisymmetric magnetic field, Begelman (1998) suggested that the problem can be alleviated 
if a current-driven (CD) kink instability destroys the concentric field structure in the nebula. 
In the axisymmetric case, magnetic loops in the expanding fiow are forced to expand and 
perform work against the magnetic tension. The CD kink instability allows the loops to 
come apart and one expects that in 3D, the mean field strength is not amplified much by 
expansion of the fiow and the hoop stress would not necessarily pinch the flow as much as 
would otherwise be supposed. In this case a just upstream of the termination shock might 
not need to be so unreasonably small as was found in axisymmetric models. 

This idea can be checked only by 3D simulations of plasma flow within the nebula. As a 
first step we simulate the 3D evolution of the simple cylindrical model of PWNe developed by 
Begelman & Li (1992). This model describes a quasi-static cylindrical configuration with a 
purely toroidal magnetic field. The plasma within the cylinder is relativistically hot and the 
hoop stress is balanced by the thermal pressure. The cylinder is confined on the outside by a 
nonmagnetized gas. The hnear analysis shows (Begelman 1998) that such a configuration is 
unstable with respect to the CD kink instability. Here we simulate the nonlinear evolution of 
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this system in order to see whether stabihzing boundary conditions suppress development of 
the instabihty inside the plasma volume via 3D relativistic magnetohydro dynamic (RMHD) 
simulations of the CD kink instability in a relativistically hot plasma column containing 
a toroidal magnetic field. This paper is organized as follows: We describe the numerical 
method and setup used for our simulations in §2, present our results in §3, and discuss the 
astrophysical implications in §4. 

2. Numerical Setup 

To study the evolution of the CD kink instability in a hot plasma column containing 
a toroidal magnetic field, we use the 3D GRMHD code "RAISHIN" in three-dimensional 
Cartesian geometry. RAISHIN is based on a 3 + 1 formalism of the general relativistic con- 
servation laws of particle number and energy-momentum, Maxwell's equations, and Ohm's 
law with no electrical resistance (ideal MHD condition) in a curved spacetime (Mizuno et 
al. 2006) 0. In the RAISHIN code, a high-resolution shock- capturing (HRSC) scheme is 
employed. The numerical fluxes are calculated using the HLL approximate Riemann solver, 
and fiux-interpolated constrained transport (fiux-CT) is used to maintain a divergence-free 
magnetic field. The RAISHIN code performs special relativistic calculations in Minkowski 
spacetime by choosing the appropriate metric. RAISHIN has proven to be accurate to sec- 
ond order and has passed a number of one- dimensional and multidimensional numerical tests 
including highly relativistic cases and highly magnetized cases in both special and general 
relativity (Mizuno et al. 2006; Mizuno et al. 2010). Here a fifth-order weighted essentially 
non-oscillatory (WENO, Jiang & Shu 1996) reconstruction scheme is built into the RAISHIN 
code and used in the simulations in order to handle turbulent structure more finely. WENO 
schemes provide highly-accurate solutions in regions of smooth flow and non-oscillatory tran- 
sitions in the presence of discontinuous waves by combining different interpolation stencils 
of order r into a weighted average of order 2r — 1. 

We consider a hot plasma column containing a pure toroidal magnetic field with radius 
R and the height being in hydrostatic equilibrium with a hot static medium. We begin 
with the cylindrical model for PWNe proposed by Begelman & Li (1992) who calculated 



^ Constained transport schemes are used to maintain divergence-free magnetic field in the RAISHIN code. 
This scheme requires the magnetic field to be defined at the cell interfaces. On the other hand, conservative, 
high-resolution shock capturing schemes (Godonov-type scheme) for conservation laws require the variables 
to be defined at the cell center. In order to combine variables defined at these different positions, the 
magnetic field calculated at the cell interfaces is interpolated to the cell center and as a result the scheme 
becomes non-conservative even though we solve the conservation laws (Komissarov 1999). 
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the axisymmetric structure of a bubble blown by a continuously injected relativistically 
hot, magnetized plasma in a surrounding medium. According to this model, the radial gas 
pressure and toroidal magnetic field profiles in the hot plasma column are given by 

Stt 4773 ' ^ ^ 

where x — r/R, and 77 is found for any x from the equation 

V + Ix'-V'^'. (3) 

In this solution, the magnetic hoop stress is balanced by the plasma pressure. Thus, the 
plasma pressure at the axis is larger than the total (plasma and magnetic) pressure at the 
edge of the bubble. It is this discrepancy that results in the elongation of the nebula. At 
X > 1, the hot plasma column is surrounded by a hot static unmagnetized medium with 
constant pressure equal to the total pressure at the bubble edge 

Pout = — , (4) 

Vo 

where rjQ = 9/4 is the solution of Eq. (3) at a; = 1. We choose p = 1 throughout the region. In 
the simulations, we set po = lO^pc^, so that the plasma is relativistically hot, i.e. pc^ << p, 
where c is speed of light. The equation of state is that of an ideal gas with p = (T — l)pe, e 
is the specific internal energy density and the adiabatic index is set to F = 4/3. The specific 
enthalpy is /i = 1 + e/c^ -\-p/ pc^. 

The above configuration was shown to be neutrally stable with respect to axisymmetric 
perturbations but unstable to the helical modes (Begelman 1998). Therefore, we choose an 
initial small radial velocity perturbation with the form given by 

5v _^ ^ ak-,x,y . f '2nkz 



•x,y/c ^—e'^Yl ^ '^V ' 



fe=l 

where is the total number of modes, (f)k is a random phase selected from the range < 
(pk < Stt, and ak;x,y is a random direction with a^.^, + a^.y — 1. We choose 5v — 0.01c in the 
simulations. In order to evaluate the effect of different initial perturbations we consider two 
cases. In case A, we choose N — 2, (p^ — 0, {ai-^, 0'i;y) — (1; 0), and (a2;x, o,2;y) — (0, 1), and 
in case B, we choose A^ = 6 and random phases and directions. The computational domain 
is a Cartesian {x, y, z) box of size QR x QR x Lz {L^ — R) with grid resolution R — 60AL 
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and Lz — 60AL. The grid resolution is the same in all directions. The simulation grid is 
periodic along the axial (z) direction. We use a larger computational domain in the x and 
y directions to avoid boundary effects. We impose reflecting boundary conditions in the x 
and y directions to ensure that the total energy is conserved. 

3. Results 

In Figure 1 we show the evolution of the magnetized column for case A when only 
a few modes are excited. The overall evolution for the multiple modes case B is similar 
but structures are more complicated because more modes are excited. The two-dimensional 
images of the gas pressure p and the perpendicular components {By and B^) of the magnetic 
field in the xz plane at a; = and the yz plane at y = arc restricted to i? < 2 in order to 
focus on the strong interaction region. The initial small velocity perturbation excites the CD 
kink instability n — 1 mode in the x-direction and n — 2 mode in the y-direction (see Fig. 
1 &X, t — 6R/c). The radial velocity (x-component and y-component) increases with time in 
the hnear growth phase. At about t — 6R/c the CD kink instability shifts to the non-linear 
phase (see also Fig. 2). In the non- linear phase, the two modes {n — 1 and n — 2) interact 
and lead to turbulence in the hot plasma column. The turbulent motion creates complicated 
magnetic field structures in the hot plasma column. Note that the gas pressure within the 
column, which was initially high in order to balance the magnetic hoop stress, decreases 
because the hoop stress weakens. Thus, as a result of the kink instability the magnetic loops 
come apart and release the magnetic stress. 

As an indicator of the growth of the CD kink instability, the time evolution of the 
volume- averaged plasma energy Ep = ph'j'^ — p and magnetic energy Em = B'^/2 for cases 
A and B are shown in Figure 2. An initial slow evolution in the linear growth phase lasts 
up to t — 6R/c for case A {t = 7R/c for case B), and is followed by a more rapid evolution 
in a nonlinear growth phase. In the linear phase the magnetic energy gradually decreases 
with a following rapid decrease in the nonlinear phase. In the nonlinear phase, this rapid 
decrease of the magnetic energy ceases at about t — llR/c ior: case A {t — lAR/c for case 
B). The same trend was seen in our previous study of the CD kink instability of a static 
cold plasma column (Mizuno et al. 2009). While the magnetic energy declines, the plasma 
energy increases gradually in the linear and rapidly in the nonlinear phase. Here growth of 
the CD kink instability leads to radial velocity increase which contributes a kinetic energy 
component to the plasma energy in both linear and non- linear phases. At about t = llR/c 
for case A {t — lAR/c for case B), increase in the plasma energy nearly ceases and the hot 
plasma column is almost relaxed. The multiple mode case B shows a more gradual evolution 



- 7- 





Fig. 1. — Two dimensional images of (case A): the gas pressure in the xz plane at y = for 
(a) t = 6R/c, (c) llR/c, the gas pressure in the yz plane at a; = for (6) t = 6R/c, (d) 
llR/c, the magnetic field By in the xz plane at y = for (e) t = 3R/c, (g) 6R/c, («) 9R/c, 
[k) 12R/c, and in the yz plane at a; = for (/) t = 3R/c, (h) 6R/c, (j) 9R/c, (/) 12i?/c. 
Arrows indicate the velocity in each plane. 
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(a) 

1.205 




Fig. 2. — Time evolution of the (a) volume-averaged plasma energy Ep, (b) magnetic energy 
Em, (c) total energy Et normalized by initial total energy {Effi), and (d) entropy for case A 
{solid lines) and case B {dashed lines). 



in the nonlinear phase and later relaxation than the two mode case A. Multiple modes lead to 
a more gradual interaction, slower development of turbulent structure, and later relaxation 
of the hot plasma column. The volume-averaged total energy Et = Ep + Em shown in Fig. 
2c is conserved at the 0.01% level in both simulations. Figure 2d shows the time evolution 
of the volume-averaged entropy S = \n{p/p'^). The entropy gradually increases in the linear 
phase and rapidly increases in the nonlinear phase. 

Our results show that the magnetic field dissipates in the nonlinear phase. This occurs 
because strong gradients in the magnetic field are formed. In our simulations, dissipation 
is numerical; in the real world, one expects the same effect because the kink instability of 
a purely toroidal magnetic field develops at very small scales and as a result dissipation is 
inevitable. The volume averaged total energy, Et = Ep + Em, of the system (Fig. 2c) shows 
that variation of the total energy is at the 0.01% level and is much less than the variation 
in the plasma energy. This indicates that our numerical dissipation conserves energy with 
sufficient accuracy. 

Figure 3 shows time evolution of the magnetization parameter a = B"^ / {ph) which 
averaged only over the magnetized hot plasma column but not over the whole volume for 
case A and B. Initially the volume-averaged magnetization cr = 0.3 in the hot plasma column. 
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Fig. 3. — Time evolution of the volume- averaged magnetization parameter a in the hot 
plasma column (-R < 1) for case A {solid line) and case B {dashed line). 

In linear growth phase, a gradually decreases in both of cases. After t = 6R/c for case A 
{t = 7R/c for case B) a rapidly decreases because the magnetic field dissipates. 

In order to study the relaxation of the hot plasma column in more detail, we show the 
time evolution of radial profiles of toroidal {(p direction) and axial {z direction) averaged 
quantities in Figures 4 and 5. Initially the hot plasma column contains toroidal magnetic 
field in hydrostatic equilibrium. In the linear phase of the CD kink instability, the radial 
and axial components of the magnetic field grow while the toroidal magnetic field and gas 
pressure decline gradually beginning near the axis. When the CD kink instability enters the 
nonlinear phase (at t = 6R/c for case A and t = 7R/c for case B), the toroidal magnetic field 
and gas pressure decrease rapidly, and the radial and axial components of the magnetic field 
increase throughout the hot plasma column. At the end of the nonlinear phase (at t = llR/c 
for case A and t = lAR/c for case B), all magnetic components become comparable and the 
field becomes totally chaotic. In a saturation phase, the magnetized column begins slow 
expansion so that eventually the column should merge with the surrounding gas. For the 
different initial perturbation profiles, the evolutionary time scale is different but the physical 
behavior is similar. Therefore, final relaxation of the hot plasma column is independent of 
the initial perturbation profile. 




4. Summary and Discussion 

We have investigated the development of the CD kink instability of a hydrostatic hot 
plasma column containing a toroidal magnetic field as a model for pulsar wind nebulae. The 
CD kink instability is excited by a small initial velocity perturbation and turbulent structure 
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Fig. 4. — Time evolution of the radial profile of the toroidal and axial- averaged radial (Br), 
toroidal (B^), and axial (Bz) magnetic field components and the gas pressure at t = OR/c 
(solid), 5R/c (dotted), 8R/c (dashed), llR/c (dash- dotted), and 15R/c (dash- double- dotted) 
for case A. 
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Fig. 5. — Time evolution of the radial profile of the toroidal and axial- averaged radial (Br), 
toroidal (-B^), and axial (Bz) magnetic field components and the gas pressure at t = OR/c 
(solid), 7R/c (dotted), llR/c (dashed), lAR/c (dash- dotted), and 18R/c (dash- double- dotted) 
for case B. 
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develops inside the hot plasma column. At the end of nonlinear evolution of the CD kink 

instability, the hot plasma column relaxes with a slow radial expansion. The magnetization a 
decreases from an initial value of 0.3 to a final value of 0.01. We find that for different initial 
perturbation profiles the time scale is a bit different but the physical behavior is the same. 
Therefore the relaxation of a hot plasma column is independent of the initial perturbation 
profile. 

Our simulations confirm the scenario envisaged by Begelman (1998). Toroidal magnetic 
loops come apart, the hoop stress declines and the pressure difference across the nebula 
is washed out. In our simulations, the ratio of the gas pressure on the axis to the total 
(magnetic+gas) pressure at the plasma column boundary decreases from 2.5 to 1.5 during 
the linear phase, while magnetic dissipation is still small. In the nonlinear phase, the mag- 
netic field dissipates and the gas pressure excess near the axis disappears. For this reason, 
elongation of a pulsar wind nebula cannot be correctly estimated by axisymmetrical models, 
because axisymmetric models retain a concentric toroidal magnetic field geometry. 

Radiation from the Crab nebula is highly polarized along the axis of the nebula, which 
is indicative of a toroidal magnetic field. We see that even though the instability eventually 
destroys the toroidal structure, the magnetic field becomes completely chaotic only at the 
end of the nonlinear stage of development. Therefore, the toroidal magnetic field should 
dominate in the central parts of the nebula that are filled by recently injected plasma. 

Our simple model does not allow us to determine the value of a in the pulsar wind; for 
this purpose one has to perform fully 3D simulations of the expanded nebula taking into 
account the continuoiis injection of plasma. Here we have demonstrated that 3D effects 
are crucially important to a determination of the structure of pulsar wind nebulae and 
that previous dynamical arguments concluding that a must be extraordinarily small can be 
abandoned. 
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